
*******************
global E_list="-3100(25)-125"
global band=100
local lhs "cog"
local rhs "hours"
local robust "simple"
local K=$Kmain
local controls "${CONTROLlist_`lhs'}"
local subsample "$subsample"
local type "hom"


qui {
	local countE=0
	forvalues E = $E_list {
		local countE=`countE'+1
	}
	local rowcounter = 0
	matrix actual_estimates_mat = J(`countE',11,-9)
	use "$datapath/003_cluster_`subsample'_`lhs'.dta", clear
	gen bunch=(hours==0)
	noisily sum bunch
	local m0=r(mean)
	noisily ccn_ESens_het het_uniform 1 `rhs' `lhs' `robust' "`controls'" het_tobit -100
	matrix actual_estimates_mat[1,1] = r(E_uniform)
	
	qui ccn_ESens_het het_tobit 1 `rhs' `lhs' `robust' "`controls'" het_tobit -100
	matrix actual_estimates_mat[1,2] = r(E_het_tobit)
					
	qui ccn_ESens_het symmetric 1 `rhs' `lhs' `robust' "`controls'" het_tobit -100
	matrix actual_estimates_mat[1,3] = r(E_symmetric)
					
	noisily ccn_ESens_het nohighpeaks 1 `rhs' `lhs' `robust' "`controls'" het_tobit -100
	matrix actual_estimates_mat[1,4] = r(E_NoHighPeaks)

	noisily ccn_ESens_het bilog 1 `rhs' `lhs' `robust' "`controls'" het_tobit -100
	matrix actual_estimates_mat[1,5] = r(E_bilog_lb)
	matrix actual_estimates_mat[1,6] = r(E_bilog_ub)						

	forvalues E = $E_list {
		noisily display "LHS=`lhs', RHS=`rhs', SE=`robust', K=`K', $S_DATE, $S_TIME, E=`E'"
		local rowcounter = `rowcounter'+1
		use "$datapath/003_cluster_`subsample'_`lhs'.dta", clear
		qui tab X_`K', generate(dumX_`K'_)
		gen A=afqt_mom
		gen AL=A*`rhs'
		qui ccn_ESens_het Ecounterfactual `K' `rhs' `lhs' `robust' "`controls'" het_tobit `E'
		matrix actual_estimates_mat[`rowcounter',7] = r(BETA)*100
		matrix actual_estimates_mat[`rowcounter',8] = r(BETA_se)*100		
		matrix actual_estimates_mat[`rowcounter',9] = r(BETA_L)*1000
		matrix actual_estimates_mat[`rowcounter',10] = r(DELTA)*100			
		matrix actual_estimates_mat[`rowcounter',11] = `E'
	}
	mat colnames actual_estimates_mat = E_uniform E_het_tobit E_symmetric E_NoHighPeaks E_bilog_lb E_bilog_ub b_E se_E bL_E d_E E
	clear
	qui svmat actual_estimates_mat, names(col)
	gen LHS="`lhs'"
	gen RHS="`rhs'"
	gen Kcluster=`K'
	gen Robust="`robust'"
}
save "$datapath/000_`type'_ESens.dta", replace
replace E_uniform = . if E_uniform==-9
replace E_het_tobit = . if E_het_tobit==-9
replace E_symmetric = . if E_symmetric==-9
replace E_NoHighPeaks = . if E_NoHighPeaks==-9
replace E_bilog_lb = . if E_bilog_lb==-9
replace E_bilog_ub = . if E_bilog_ub==-9
sum E_uniform
local E_uniform = r(mean)
sum E_het_tobit
local E_het_tobit = r(mean)
sum E_symmetric
local E_symmetric = r(mean)
sum E_NoHighPeaks
local E_NoHighPeaks = r(mean)
sum E_bilog_lb
local E_bilog_lb = r(mean)
sum E_bilog_ub
local E_bilog_ub = r(mean)

keep b_E se_E bL_E d_E E*
tempfile temp
save `temp', replace

use "$datapath/001_nlsy_child_mother_cog.dta", clear
drop if hours==2080
keep if hours<=2400
local mainRHS = "hours"
local skill="cog"
local robust="simple"
gen X_1=1

qui ccn_fit tobit  `mainRHS' `skill' 100 `robust' ""
qui ccn_fit heterotobit  `mainRHS' `skill' 100 `robust' ""
qui ccn_fit heterounif  `mainRHS' `skill' 100 `robust' ""
qui ccn_fit heterosymm  `mainRHS' `skill' 100 `robust' ""
qui ccn_fit heterosymmetry  `mainRHS' `skill' 100 `robust' ""
qui ccn_fit nohighpeaks  `mainRHS' `skill' 100 `robust' ""
qui ccn_fit bilog  `mainRHS' `skill' 100 `robust' ""


use "$datapath/002_Fit.dta", clear
tempfile temp2
save `temp2', replace

use "$datapath/001_nlsy_child_mother_cog.dta", clear
drop if hours==2080
keep if hours<=2400
local rhs = "hours"
sort hours
cumul hours, generate(CDF_hours)
label variable CDF_hours "All Sample"			
append using `temp'
append using `temp2'
cumul muH_HeteroUnif, generate(CDFmuH_HeteroUnif)
replace CDFmuH_HeteroUnif=0 if CDFmuH_HeteroUnif==.
reg muH_HeteroUnif CDFmuH_HeteroUnif if CDFmuH_HeteroUnif>=0.01
predict muhat_HeteroUnif

gen CDFmuH_HeteroTobit=.
replace CDFmuH_HeteroTobit=q/100 if q>0

/*
cumul muH_HeteroTobit, generate(CDFmuH_HeteroTobit)
replace CDFmuH_HeteroTobit=0 if CDFmuH_HeteroTobit==.
reg muH_HeteroTobit CDFmuH_HeteroTobit if CDFmuH_HeteroTobit>=0.01
predict muhat_HeteroTobit
*/

cumul muH_NoHighPeaks, generate(CDFmuH_NoHighPeaks)
replace CDFmuH_NoHighPeaks=0 if CDFmuH_NoHighPeaks==.
reg muH_NoHighPeaks CDFmuH_NoHighPeaks if CDFmuH_NoHighPeaks>=0.01
predict muhat_NoHighPeaks

replace CDFmuH_bilog_lb=. if CDFmuH_bilog_lb<0
replace CDFmuH_bilog_ub=. if CDFmuH_bilog_ub<0
gen muH_bilog = q

replace E_uniform=`E_uniform'
replace E_het_tobit=`E_het_tobit'
replace E_symmetric=`E_symmetric'

local E_uniformT=`E_uniform'+20
local E_het_tobitT=`E_het_tobit'+20
local E_symmetricT=`E_symmetric'+20
local E_NoHighPeaksT=`E_NoHighPeaks'+20

twoway (line b_E E if E>-3001, lpattern(solid) lwidth(thick) lcolor(black)) || ///
		, scheme(s1color) ylabel(-.1(.02)0, axis(1)) xlabel(-3000(1000)0) ytitle("Estimate of {&beta}") xtitle("E [ L* | L = 0 ]") ///
		xline(`E_uniform', lpattern(dash) lcolor(blue)) xline(`E_het_tobit', lpattern(dot) lcolor(green)) xline(`E_symmetric', lpattern(solid) lcolor(black)) ///
		 text(-.105 `E_uniformT' "E{sub:U}", color(blue)) text(-.105 `E_het_tobitT' "E{sub:N}", color(green)) text(-.105 `E_symmetricT' "E{sub:S}") 
graph export "$figpath/b_ESens_hom.pdf", replace



twoway (line CDF_hours hours, lcolor(black) lwidth(thick) sort)  ///
		|| (line CDFmuH_NoHighPeaks muhat_NoHighPeaks if muhat_NoHighPeaks<=0, lpattern(longdash) lcolor(red) sort) ///
		, legend(off) scheme(s1mono) ylabel(0(.2)1) ///
		xlabel(-1000(1000)2000) ytitle(CDF) xtitle("           E [ L* | L = 0 ]                                                                                             L=L*  ") xline(`E_uniform', lpattern(dash) lcolor(blue)) ///
		xline(`E_het_tobit', lpattern(dot) lcolor(green)) xline(`E_symmetric', lpattern(solid)) ///
		xline(`E_NoHighPeaks', lpattern(longdash) lcolor(red))	///
		text(-.055 `E_uniformT' "E{sub:U}", color(blue)) text(-.055 `E_het_tobitT' "E{sub:N}", color(green)) text(-.055 `E_symmetricT' "E{sub:S}") ///
		text(-.055 `E_NoHighPeaksT' "E{sub:0}", color(red))
graph export "$figpath/CDF_ESens_hom.pdf", replace



